Reaction-controlled diffusion: Monte Carlo simulations 
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We study the coupled two-species non-equilibrium reaction-controlled diffusion model introduced 
by Trimper et al. [Phys. Rev. E 62, 6071 (2000)] by means of detailed Monte Carlo simulations 
in one and two dimensions. Particles of type A may independently hop to an adjacent lattice site 
provided it is occupied by at least one B particle. The B particle species undergoes diffusion-limited 
reactions. In an active state with nonzero, essentially homogeneous B particle saturation density, 
the A species displays normal diffusion. In an inactive, absorbing phase with exponentially decaying 
B density, the A particles become localized. In situations with algebraic decay ps(f) ~ t~ aB , as 
occuring either at a non-equilibrium continuous phase transition separating active and absorbing 
states, or in a power-law inactive phase, the A particles propagate subdiffusively with mean-square 
displacement (x(t) A ) ~ t 1 ~ aA . We find that within the accuracy of our simulation data, q.a ~ as as 
predicted by a simple mean-field approach. This remains true even in the presence of strong spatio- 
temporal fluctuations of the B density. However, in contrast with the mean-field results, our data 
yield a distinctly non-Gaussian A particle displacement distribution nA(x,t) that obeys dynamic 
scaling and looks remarkably similar for the different processes investigated here. Fluctuations of 
effective diffusion rates cause a marked enhancement of UA(x,t) at low displacements |a?[, indicating 
a considerable fraction of practically localized A particles, as well as at large traversed distances. 

PACS numbers: 05.40.-a, 05.40.Fb, 64.60.Ht, 82.20.-w 



I. INTRODUCTION 



The goal of statistical mechanics is to understand the 
relationship between microscopic and macroscopic dy- 
namics in systems consisting of a large number of de- 
grees of freedom. One classical success of the equilibrium 
formalism is the prediction of universal phase transition 
behavior: independent of the microscopic details of their 
interactions, systems with identical overall features, gov- 
erned by their symmetries, spatial dimension d, and per- 
haps large-scale interaction properties, display very sim- 
ilar phase diagrams. Moreover, their critical points are 
characterized by the same small set of independent scal- 
ing exponents. Thus physical systems with very compli- 
cated interactions can often be adequately described by 
considerably simplified models, which in turn form the 
basis of simulation studies and numerical analysis. 

Here we investigate a simple coupled reaction-diffusion 
system, which however leads to remarkably rich features. 
More specifically, the spatio-temporal fractal structures 
emerging at a non- equilibrium critical point of a reacting 
species B impose non-trivial scaling behavior onto the 
propagation of passive random walkers A, whose propa- 
gation is however limited to sites occupied by at least one 
B particle. One may envision this system to model virus 
(represented by the A particles) propagation in a carrier 
B population that is set close to its extinction threshold; 
the virus remains dormant when there are no B organ- 
isms present. Below, we shall encounter and character- 
ize the ensuing scaling laws by means of Monte Carlo 
simulations, and compare our numerical results with the 



predictions of a mean-field approximation. 

In non-equilibrium systems, the detailed balance con- 
ditions are violated; i.e., the probability of at least one 
closed loop of transitions between microscopic configura- 
tions depends upon the direction the loop is traversed. 
This is the case even in stationary states in open sys- 
tems, through which a steady particle or energy current 
from the outside is maintained. Outside physics non- 
equilibrium models may describe, for example, popula- 
tion dynamics, chemical catalysis, and financial markets. 
Yet reassuringly, universal behavior has also been found 
to persist for non-equilibrium models that display phase 
transitions between different stationary states. 

Prominent examples are continuous transitions be- 
tween active and inactive/ absorbing states in diffusion- 
limited 'chemical' reactions [1J. The class of models we 
will be studying involves competing annihilation and off- 
spring reactions of a single species B, performing unbi- 
ased random walks on a d-dimensional hypercubic lattice: 
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with integers m > 1, n > 1. For large branching rate cr, 
the system is in an active state with a non-zero and es- 
sentially homogeneous particle density. In contrast, when 
the decay processes with rates A and n dominate, the B 
particle density reaches zero, and the dynamics ceases en- 
tirely in this inactive, absorbing state. By appropriately 
tuning the reaction rates a continuous phase transition 
between these two stationary states can be observed Q . 
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Generically, such transitions fall into the directed per- 
colation (DP) universality class [3j with upper critical 
dimension d c = 4. The standard example is represented 
by the Gribov process B — + 0, B =^ 2B. Equivalcntly, 
one may use the scheme Q with m — 1 and n — 2. Di- 
rected percolation was initally devised to characterize the 
transition from finite to infinite-sized clusters in directed 
media (such as a porous rock in a gravitational field) £| . 
Other applications include certain models of catalytic re- 
actions, interface growth, turbulence, and the spread of 
epidemics (5j. Experimental evidence for DP critical be- 
havior was recently observed in spatio-temporal intermit- 
tency in ferrofluidic spikes [||. 

For active to absorbing state phase transitions in 
single-species reactions that include first-order processes, 
in the absence of memory effects and quenched disorder, 
the so-called parity conserving (PC) universality class ap- 
pears to represent the only scenario for non-DP critical 
behavior In the above reaction scheme (JJJ, PC scal- 
ing is observed for branching and annihilating random 
walks with A = 0, n — 2, and even offspring number to. 
In that situation, reactions either create or annihilate an 
even number of particles. Thus the number of B parti- 
cles remains either even or odd throughout the system's 
temporal evolution. Indeed, the distinct non-trivial scal- 
ing exponents of the PC universality class, albeit limited 
essentially to d = 1, can be attributed to this special 
symmetry of local particle number parity conservation. 
Moreover, for d < d' c w 4/3, fluctuations cause the emer- 
gence of a power-law inactive phase, characterized by the 
algebraic decay laws of diffusion-limited pair annihilation 
25^0 (A = a = 0, n = 2) 0. 

When A > or m is odd (for n = 2), however, parity 
conservation is destroyed. The case A > immediately 
yields a transition in the DP universality class with d c — 
4. Yet for odd to, even if A = initially, fluctuations 
generate sufficiently strong decay processes B — ► in 
d < 2 dimensions to produce a transition to an inactive 
phase with DP critical exponents. For A = and d > 2, 
one encounters only an active phase for any a > 0, as 
predicted by the mean- field rate equation Q- 

At the non-equilibrium continuous phase transition, 
the reacting particles form spatio-temporal fractal struc- 
tures characterized by algebraic decay of the correlation 
functions (one example is depicted in Fig.QJ. In Ref. @ it 
was suggested to employ these dynamic fractals of react- 
ing agents B as backbones for nearest-neighbor hopping 
processes of another, otherwise passive particle species 
A. The A particles are then allowed to move to a nearest 
neighbor site only if that site is occupied by at least one B 
particle. In an active state of the B system, with a largely 
homogeneous particle distribution, the A particles follow 
Fick's normal diffusive propagation law, (x(t) A ) =2Dt, 
with diffusion constant D ~ Oq/to, where do denotes 
the lattice constant of the hypercubic lattice, and Tq 
the microscopic hopping rate. On the contrary, in a DP- 
type inactive phase with an exponentially decaying B 
particle density, the A species will become localized, i.e., 
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FIG. 1: Space-time plot for the one-dimensional reactions 
B 3B, 2B — > at the active/absorbing critical point (PC 
universality class). 

(x(t)\) — > const as t — ► oo. In precisely this sense the 
inactive to active state transition of the B system thus 
induces a localization transition for the A particles. At 
the transition itself, as well as in the PC-inactive phase, 
the B density decays algebraically, 

PB(t)~t- aB , (2) 

with < otB < 1- Correspondingly, the A species propa- 
gates subdiffusively, 

(m a A)~*- aA > (3) 

where again < a a _ 1- In fact, a simple mean-field 
approach suggests a a = O-b @- Our goal here is to fur- 
ther elucidate the scaling laws for the ensuing anomalous 
A particle diffusion through Monte Carlo simulations in 
one and two dimensions. We shall also numerically de- 
termine the full time-dependent probability distribution 
for the A species displacements and compare it with the 
Gaussian distribution predicted by mean-field theory . 

Before we proceed, we note that our model is related 
to, but quite distinct from studies of anomalous diffusion 
on static fractal structures. These have been used to de- 
scribe a variety of physical phenomena, such as percola- 
tion through porous or fractured media and electron-hole 
recombination in amorphous semiconductors |f| . Anoma- 
lous diffusion may also ensue for particle transport in a 
random medium with quenched disorder, provided the 
obstacles are sufficiently strong to effectively reduce the 
number of diffusive paths on large length scales [ljj. We 
emphasize again that in the system under investigation 
here the fractal structure evolves dynamically, which pro- 
vides an alternative mechanism for subdiffusive propaga- 
tion: When Eq. (J2J holds, the number of available paths 
decreases with time. Yet notice that A particles that 
have become localized on an isolated B cluster for some 
time may still become linked to a large connected region 
of available sites later. 

In the subsequent Sec. [H] we briefly review the theoret- 
ical considerations of Ref. [8( , and list the central results 
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TABLE I: Critical exponents for the parity conserving (PC) and directed percolation (DP) universality classes of active to 
absorbing phase transitions, as determined from Monte Carlo simulations in one and two dimensions [1]. Here, r denotes the 
deviation of a relevant control parameter from the critical point. For directed percolation, the first-order results from the e 
expansion near d c = 4 dimensions are given as well. 



Critical exponents 



PC, d 



DP, d = 1 



DP, d - 



DP, d - 



e~lrr 

t c ~ £ z ~ \r 
pc{t) ~ r 



/3 » 0.92 
!/ w 1.84 
z » 1.75 
a » 0.285 



/3 « 0.276 
i/ « 1.097 
z » 1.581 
a w 0.159 



/3 « 0.584 
i/ « 0.734 
z » 1.764 
a w 0.451 



/3 = 1 - e/6 + 0(e 2 ) 
v = l/2 + e/16 + 0(e 2 ) 
z = 2 - e/12 + 0(e 2 ) 
a= l-e/4 + 0(e 2 ) 



from the mean-field approach for the A species propa- 
gation on the dynamic B fractal. In Sec. II I II we give 
an overview of the Monte Carlo simulation methods em- 
ployed in our study. Next, in Sec. IIVI we first present 
our results for the anomalous A diffusion as induced by 
pure annihilation kinetics of the B species (a — 0). Sec- 
tion E] is devoted to the central issue in our investiga- 
tion, namely the subdiffusive behavior of the A particles 
at the localization transition caused by the active to in- 
active/absorbing phase transition of the reacting agents 
B. We summarize our results in Sec. IVII add concluding 
remarks, and point out a few open problems. 



II. THEORETICAL PREDICTIONS 

For the specific combination of B particle reactions 
studied here, the asymptotic scaling behavior is well 
understood. In particular, we shall consider both sys- 
tems with pure annihilation kinetics (n B — > 0) , and 
models with competing annihilation and offspring re- 
actions, as listed in |JI)|. in the vicinity of their non- 
equilibrium phase transition from active to inactive, ab- 
sorbing states. These phase transitions either fall into 
the directed percolation (DP) or parity conserving (PC) 
universality classes, both of which have been previously 
investigated in detail 0. Specifically, the exponents a B 
characterizing the long-time decay of the particle density 
are well-established through simulations (see Tab.QJ. 

Let us first consider pure annihilation kinetics with re- 
action rate p. The corresponding mean-field rate equa- 
tion for the B particle density reads 



dt PB{t) = -np,p B (t) 7 
For n > 1 this yields 



PBit) 



Pb(0) 



(l + f/O 1 ^ 



p B (oy 



n(n — 1) p 



(4) 



, (5) 



whence at sufficiently large times p B (t) ~ (/ii)~ 1 / ( -™~ 1 \ 
independent of the initial density p B (0). For n = 1, i.e., 
spontaneous death with rate A, one naturally finds expo- 
nential decay, 



We expect these results to be valid (at least quali- 
tatively) in dimensions above an upper critical dimen- 
sion d c , which can be determined through straightfor- 
ward dimensional analysis: since the B particles un- 
dergo ordinary diffusion, we expect [t] = [x] 2 , and of 
course in d dimensions, [p B ] = [x]~ d - Then by Eq. 



[/<■] 



ld(ra-l)-2 



Thus (i is dimensionless (marginal in 



the RG sense) at d c (n) = 2/(n — 1). The mean-field 
description J2J| fails in lower dimensions where there is 
a non-negligible probability that a particle will retrace 
part of its trajectory. For n — 2, 3 anti-correlations be- 
tween surviving particles are induced in dimensions d < 2 
and d = 1, respectively, since many of the nearby par- 
ticles along a specific agent's trajectory are annihilated 
in the first trace. The annihilation processes then be- 
come diffusion- rather than reaction-limited. For the pair 
process, the particles need to traverse a distance £(t) ~ 
{Dt) 1 / 2 before they can meet and annihilate, where D de- 
notes the B particle diffusion constant. Hence the parti- 
cle density should scale as n{t) ~ l{ty d - (Dt)~ d / 2 . 
Indeed, a renormalization group analysis predicts for 
2B — > (as well as for pair coagulation 2B — > B) 



PB (t) ~ (Dt)- d / 2 
p B (t) - {Dt)- 1 In Dt 



for d < 2 , 
at 4(2) = 2 , 



(7) 
(8) 



in agreement with exact solutions in d = 1. Thus parti- 
cles survive considerably longer than Eq. JSJ would sug- 
gest. For triplet annihilation (3£> — > 0), in one dimension 
there remain mere logarithmic corrections to the mean- 
field result [13, 



Ps{t) 



InDt 
Dt 



1/2 



at d c {3) = 1 



(9) 



p B (t) = p B (0) e 



-At 



(6) 



The higher-order (n > 4) annihilation processes should 
all be aptly described by the mean-field power laws 10). 

Exact results for the critical behavior of the DP and 
PC universality classes cannot be derived analytically, 
but the scaling exponents have been measured quite ac- 
curately by means of computer simulations see Tab.[I] 
The universal properties of directed percolation can be 
represented through Reggeon field theory which al- 
lows a systematic perturbational calculation of the crit- 
ical exponents in an e expansion near its upper critical 
dimension d c — 4. The one-loop fluctuation corrections 
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to their mean-field values are listed in Tab. [I] as well. 
At least for DP, the scaling relation (3 — z v a holds. A 
similarly reliable analytic computation of the PC critical 
exponents has as yet not been achieved, owing to the ab- 
sence of a corresponding mean-field theory (see Ref. |?J 
for further details). 

In Ref. 0, the reaction-controlled diffusion model was 
defined as follows. Otherwise passive agents A perform 
independent random walks to those sites that are occu- 
pied by at least one B particle; the B species is subject 
to diffusion-limited reactions of the above type. In our 
simulations, we have assumed the A hopping rate Tq 1 
to be independent of the number of B particles at ad- 
jacent sites. This contrasts the model of Ref. @, where 
the A hopping probability to a given site was taken to 
be proportional to the number of B particles on that 
site. Yet as we are mostly interested in the asymptotic 
behavior at low B densities, this distinction should be 
largely irrelevant. Moreover, in the majority of systems 
studied here the B pair annihilation reaction was set to 
occur with probability 1, which eliminates multiple site 
occupations. 

To be specific, consider B particles undergoing the Gri- 
bov reactions B — > 0, B ^ 2B, with an ensuing critical 
point in the DP universality class. The effective theory 
near the phase transition then becomes equivalent to a 
Langevin equation for a fluctuating field b(x, t) 0,0,^3 : 

d t b = D (V 2 - r) b~2pb 2 + r 1 . (10) 

With (r](x,t)) = 0, the ensemble average of b over noise 
realizations yields the mean B particle density, {b(x, t)) = 
PB(t)- For the correlator of the stochastic noise one finds 

(rj(x, t) r)(x', t'))=2a b(x, t) 5{x - x') 5(t - t') , (11) 

which is to be understood as the prescription to always 
factor in the local particle density when noise averages 
are taken. According to Eq. (|11|) . all fluctuations vanish 
in the absorbing state, as they should. As before, A, 
a, and /i represent the B particle decay, branching, and 
coagulation rates, respectively. The control parameter r 
denotes the deviation from the critical point, e.g., r = 
(A — cr)/D in the mean-field approximation. 

With the model definition in Ref. Q , the effective diffu- 
sivity of the agents A becomes proportional to the local B 
density. In fact, starting from the classical master equa- 
tion, one can derive the following continuum stochastic 
equation of motion for a coarse-grained field a{x 1 1) that 
describes the A species Q 

d t a = D{V 2 a)b~ Da{V 2 b) + C , (12) 
with noise correlations 

(((x,t)((x / ,t')) = J (13) 
«(x, t) r,{x', t')) = D [V 2 a(f, t)] b(x, t) S(x ~ x') 5(t - t') 
—D a(x, t) V 2 [b(x, t) S(x - x' ) S(t - t')] . 



The fluctuations of the b field thus influence the A diffu- 
sion in a non-trivial manner. 

Certainly outside the critical regime, well inside either 
the active or inactive phases, which for DP are both char- 
acterized by exponentially decaying correlations in space 
and time, one may apply a mean-field type of approxi- 
mation. To this end, we consider the B particle density 
to be spatially homogeneous, and neglect the noise cross- 
correlations in Eq. i|13fl . Upon rescaling, the equation of 
motion H12|l then reduces to a mere diffusion equation 8] 

d t n A (x,t) = Dp B {t)W 2 n A {x,t) (14) 

for the probability n A of finding a particle A at point x at 
time t, with time-dependent effective diffusivity D 
We may interpret this result as follows. In the original 
model, the hopping rate to an adjacent site is propor- 
tional to the number of B particles on that site. At suffi- 
ciently low densities that multiple B particle occupation 
of a given site can be neglected, the average B density 
represents the fraction of lattice sites available for the A 
particles to hop to. Thus we expect the A species diffu- 
sion rate to be approximately proportional to the global 
B particle density. For this assumption to be accurate, 
the A particle distribution must also be assumed to be 
at least roughly uniform. However, local fluctuations of 
the B density may induce some clustering for many A 
particles as well, though certainly to a lesser degree since 
in regions with small disjoint B clusters, the A particles 
become localized. Any deviations from the uniform effec- 
tive diffusion coefficient caused by an inhomogeneous B 
particle distribution would be diminished by this simul- 
taneous clustering. We will discuss these effects in more 
detail below as we measure deviations from the mean- 
field behavior. 

Under the above mean-field assumption, we may em- 
ploy the diffusion equation (|14|) to determine how the 
probability distribution n A (x, t) evolves in time for a sys- 
tem of independent A particles that all start initially at 
a particular location x = 0, i.e., n A (x, 0) = S(x). Even 
with a time-dependent diffusion coefficient, Eq. (|14Jl is 
readily solved via spatial Fourier transformation, result- 
ing in a Gaussian: 

rU(f '^i47^^ eXP ("^) ■ (15) 

But the expression Dt in Fick's law for standard diffusion 
becomes replaced with an integral over the evolving B 
density, 

D'{t) = D f PB (t')dt' . (16) 
Jo 

Naturally, the odd moments of the distribution (|15fl van- 
ish, while (x(t) k A ) = (\x(t) A \ k ) > for k even. We then 
compute 

(\x(t) A \ k ) = f \x\ k n A (x,t)d d x 
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T(d/2) { 2 J ■ [ '> 

For k = 2 in particular, this reduces to 

(x(t) 2 A ) = 2dD'(t) . (18) 

For a constant B particle density, e.g., the saturation 
value p s in an active phase, we recover ordinary diffusion 
with effective diffusivity D — Dp s . In an inactive phase 
with exponential density decay I©, we find instead 

(^) = ^f^(l-e-^) • (19) 

Thus, asymptotically the A particles become localized, 
with (x(t) 2 A ) -»■ 2dDp B (0)/X in the limit t -> oo. On 
the other hand, if the B density decays algebraically, 
see Eq. then for < a B < 1 our mean-field solu- 
tion predicts subdiffusive propagation with a a = olb-, 
whereas asymptotically 

(f(i)i) ~ D InDt (20) 

if as — 1. Finally, for pair annihilation processes at the 
critical dimension d c = 2, governed by Eq. (jSJ), one finds 

(x(t) 2 A )~D(lnDt) 2 , (21) 

and similarly we obtain for triplet annihilation in one 
dimension, see Eq. J§J|, 

(x(t) 2 A ) ~ D[Dt InDt) 1 / 2 . (22) 

III. MONTE CARLO SIMULATION METHODS 

Our goal was to employ Monte Carlo simulations to 
compare a a with a B , as well as to determine the dis- 
placement probability distribution n A (x, t) for the A par- 
ticles, and look for deviations from the Gaussian distri- 
bution (|15|l predicted by the mean-field approach. The 
simulations discussed here were executed on a cubic lat- 
tice with periodic boundary conditions in each spatial di- 
rection. In one dimension the lattice contained between 
10 4 and 10 5 sites, and the two-dimensional lattice ranged 
in size from 100 x 100 to 800 x 800. In each simulation 
the system was initialized by putting one B particle at 
each site (initial density p B (0) = 1 0]), and by ran- 
domly placing throughout the lattice a fixed number of 
A particles with no site exclusion. For all of the data 
given below (except where otherwise noted) , the density 
of A particles in the lattice was fixed at pa = 0.5. Each 
time step involved a complete update of the A species, 
followed by a complete update of the B particles. The 
simulation was terminated either when the number of B 
particles reached a certain lower limit (usually 0.1% of 
the number of lattice sites), or after a fixed number of 
time steps (typically « 10 6 ). 



With the exception of some particular runs in which 
the B particles were non- interacting (i.e., subject only to 
the decay B — > 0) and could thus be updated serially, we 
proceeded as follows. Given N particles at the beginning 
of the time step, N random B particles on the lattice 
were chosen to be updated. In a given time step, some B 
particles might then be addressed more than once while 
others not at all. This is appropriate even though the 
number of B particles is changing in time, because the 
net loss of B particles per time step becomes less than one 
on time scales short compared to the simulation length. 
The B particles were in general subjected to the reac- 
tions listed in and a single update proceeded in that 
sequence. To implement the process B — > 0, the B par- 
ticle was deleted with probability A. Next, the B parti- 
cle underwent an offspring reaction B — > (m + 1) B with 
probability a. In the simulations discussed here, we chose 
m = 1, 2, or 4. The offspring particle(s) were placed on 
the parent particle's nearest neighboring sites such that 
no offspring were placed on identical places, and were 
then subject to the annihilation reaction nB — > (with 
n = 2, 3, or 4) if applicable. If the B particle being 
updated did not undergo an offspring reaction, it sub- 
sequently hopped to a nearest neighbor site with some 
probability and was then subject to annihilation, which 
required all n B particles to be located on the same site. 

Initially, the A particles were also updated via Monte 
Carlo. To this end, a random direction was chosen, and 
the A particle was moved one step in that direction pro- 
vided at least one B particle was present on the desti- 
nation site. However, since the A particles were inde- 
pendent, it was later determined that they could be pro- 
cessed serially (i.e., by simply passing through the list 
of A particles) . Technically this represents a microscopi- 
cally different method, since the Monte Carlo procedure 
causes some particles to be updated more than once, and 
others not at all at a given time step. Yet we found that 
both variants produced identical macroscopic results. 

Our main interest was to determine the asymptotic 
scaling behavior of the global B particle density in the 
lattice, as well as to measure the mean-square displace- 
ment of the A species, both as functions of time t. Hence- 
forth, lengths will be measured in units of the lattice 
constant ao, and time in units of Monte Carlo steps. In 
most of our simulations we expect spatial inhomogene- 
ity in the B particle distribution. In particular, anti- 
correlations should develop in low dimensions for pure B 
annihilation kinetics. At the critical point for systems ex- 
hibiting phase transitions, at sufficicently large times the 
B species distribution should become a scale-free spatial 
fractal at length scales large compared to the lattice con- 
stant and smaller than the system size (compare Fig.^l. 
Therefore we also periodically recorded the coordinates 
of both A and B particles in order to compute probabil- 
ity distributions and correlation functions. For example, 
the B density correlation function is defined as 

C B (x,t;x',t') = (p B {x,t)p B {x',t')) - (p B ) 2 . (23) 
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By the translational and rotational invariance of the lat- 
tice, the equal-time density correlations are really a func- 
tion of la? — x'\ only. At criticality, we have 



C{x,t-x\t)=C{\x~x'\) ~ \x-x'\- 2p/u 



(24) 



whence we find at equal positions 

C(x, t; S, t 1 ) =C(\t- t'\) ~ \t - t'\~ 2f3/zu , (25) 

as a consequence of dynamic scaling. 

To compute the equal-time correlation function l|24|) 
in one dimension numerically, we fix a particular B site 
and observe the distribution of all other B particles as 
a function of the distance from it. We measure this dis- 
tribution for each B particle fixed and then average the 
resulting distributions to obtain (pB{x,t) pb(x' , i)) . In 
higher dimensions, we use the lattice directions as repre- 
sentative of the full distribution and compute C(\x — x'\) 
for pairs (x 7 x') with a common lattice coordinate in one 
direction. 

In many situations we expected the measured quanti- 
ties to be power laws as a function of time t. The simplest 
approach to computing the exponent a in a power law 
relationship p oc t~ a naturally is linear regression on In p 
vs. hit. At the continuous phase transition separating 
active and absorbing states, one expects such a power 
law dependence. However, when the system is slightly 
above or below the critical parameters, it usually behaves 
critically for some time before crossing over to super- or 
subcritical behavior (either an exponential approach to 
the final B particle density, or a power law with a dif- 
ferent exponent a'). If the system is not precisely at 
criticality (at least for the time scales simulated), then 
offcritical behavior could lead to incorrect determination 
of the critical exponent via linear regression. However, 
one can also compute a local exponent a b for a measured 
quantity p given by the expression 



a b = -log b2 [p(bt)/p(t/b)} 



(26) 



Thus at time t, supposing a power law dependence on t, 
a b is the exponent inferred from the values of p at b t and 
t/b; i.e., these two data points define a line on a log- log 
plot whose slope is — a b . 

The most time-consuming procedure in this study in- 
volved finding the critical parameters for which the sys- 
tem was at criticality. Typically, the parameters A and p, 
were fixed, and a was varied. The critical value of a was 
deduced by simultaneously increasing the lower bound 
by definitely identifying systems as subcritical, and de- 
creasing the upper bound by characterizing supercritical 
systems. In either case, the system behaved critically for 
some time (longer for a closer to the true critical value 
cr c ) before crossing over to its asymptotic behavior, and 
so critical power laws could be approximated from the 
system's intermediate scaling behavior. As noted pre- 
viously, the B species phase transition between active 
and absorbing states induces a localization transition for 



the A particles, with critical subdiffusive behavior. The 
phase transitions for both A and B particles clearly occur 
at the same value of cr c , and hence a c can be determined 
independently by measuring both the B particle density 
decay and the A species mean-square displacements as a 
function of time. We estimate our typical errors in deter- 
mining critical exponents and the subdiffusive A species 
power laws to be « ±0.01. 

We also attempted to improve the precision of our es- 
timation of <t c by a method suggested in Ref . , which 
we now briefly describe. For the measured quantity p, 
which at the critical parameter value a c depends only 
on some power of time t, i.e., p(a c ,t) ~ t~ a with some 
exponent a, one expects in the critical regime a w a c : 



p(a,t) ~ p(a c ,t) \l + ct 1/ut (a-a c ) 



(27) 



Here v t = z v describes the critical slowing down as the 
control paramater a approaches the phase transition at 
<7 C , r c ~ <r c \~ Ut ■ One first estimates a c , z/ t , and c from 
relatively short simulations at various values of a. A sim- 
ulation reaching large values of t is then performed at the 
estimated a c . One obtains an improved estimate for a c 
by replacing a and p(a, t) in Eq. (|27H with the long simu- 
lation measurements, and then finding the value of cr c for 
which p((7 c , t) is a straight line. This process may then be 
repeated to obtain the desired or computationally acces- 
sible accuracy. Note that inaccuracies in the estimates of 
c and v t result in second-order inaccuracies for a c and a 
15]. However, this method only works when a is already 
sufficiently close to a c so that this first-order approxima- 
tion is valid. In addition, statistical fluctuations in the 
data must be smoothed as much as possible through aver- 
aging over multiple runs for each a value so that c and v t 
can be somewhat accurately determined. Unfortunately, 
we found our simulations were not extensive enough in 
most cases to apply this method and obtain even more 
reliable estimates of c, vt, and a c . 



IV. ANNIHILATION KINETICS AND 
ANOMALOUS DIFFUSION 

We begin with the results of our Monte Carlo simu- 
lations for pure B species annihilation reactions. These 
serve to test the mean-field description of the ensuing 
A particle anomalous diffusion in dimensions below, at, 
and above the critical dimension d c (n) = 2/(rt — 1), and 
moreover provide a means to estimate the magnitude of 
errors to be expected in our data. In addition, the re- 
sults for the B pair annihilation model should describe 
the subcritical behavior for the reactions exhibiting ac- 
tive to absorbing transitions in the PC universality class. 



A. Spontaneous decay S — > 

We first verified that our simulations correctly repro- 
duced the n = 1 solution 10 to the mean-field equa- 
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tion J2J giving exponential B density decay. This re- 
sult should be valid in any dimension since all particles 
evolve independently. The mean-field description for the 
A particles then predicts their localization according to 
Eq. (|19fl . We ran simulations in both one (system size 
10000) and two dimensions (system size 100 x 100) using 
a decay rate A = 0.01, and indeed found excellent agree- 
ment at large times with the prediction i|19|) . Initially, 
however, the A particles moved slower than suggested 
by Eq. ljT§|, . consistent though with our different mi- 
croscopic realization of the reaction-controlled diffusion 
model: the rules adopted in the simulations only distin- 
guish between sites that are occupied or unoccupied by 
B particles, and so multiple site occupation effectively 
corresponds to lower local densities ps in the analytical 
description. Yet at low B particle densities multiple site 
occupation becomes negligible, leading to A species lo- 
calization precisely as described by Eq. i|19fl . We also 
computed the actual A particle displacement distribu- 
tion as a histogram of final net displacements. Using the 
measured average final mean-square displacement as in- 
put to Eq. I|15[l rather than estimating the coefficient D, 
we found good agreement with the predicted Gaussian 
distribution in both one and two dimensions, although 
in d = 1 the data perhaps indicate a slight excess of 
particles localized very near their initial location. 

In the results summarized above, at each Monte Carlo 
step the updated B particles hopped to a nearest neigh- 
boring site with probability 1. We also investigated the 
effect of varying this probability. Initially, the mean- 
square A displacement then grows faster in situations 
with low B particle diffusivity because the probability of 
multiple B site occupation is reduced, leaving a larger 
fraction of available sites for the A species. However, as 
the B particles are depleted the connectivity of the lattice 
decreases, and the movement of the B species quickly be- 
comes the dominant mechanism for A particle diffusion. 
We found an overall monotonic increase of final A par- 
ticle mean-square displacements as a function of the B 
diffusivity. In one dimension the (temporary) localiza- 
tion of an A particle requires only two sites unoccupied 
by B particles, compared with four sites in two dimen- 
sions. Therefore the connectivity decreases more sharply 
in d = 1 as a function of the B density. Consequently 
diminishing the B species random walk probability has a 
more pronounced effect in one dimension than in d = 2. 



B. Pair annihilation 2B -> 

Since the critical dimension for the B particle pair an- 
nihilation reaction is d c (2) = 2, the mean-field prediction 
does not apply to either the one- or two-dimensional 
simulations executed. In d = 1 we expect according to 
Eq. 10 p B (t) ~ i -1 / 2 , and using Eqs. (H§ and (dJ) the 
mean-field approach suggests (x(t)\) ~ t 1 / 2 . Figure [21 
shows the simulation results (averaged over 20 runs) for 
B pair annihilation in one dimension (system size 10000) 
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FIG. 2: The B particle density and corresponding A species 
mean-square displacement for the pair annihilation reaction 
2 B — > with fi = 0.5 in d = 1 demonstrating excellent 
agreement with the predicted asymptotic power laws. 



with annihilation probability p = 0.5. At first the B 
density decays faster than Eq. (Q predicts: initially the 
effective power law should be close to the mean-field re- 
sult ps(t) ~ t' 1 since the anti-correlations described in 
Sec. Ull develop only after some time has elapsed, where- 
upon the asymptotic decay ~ t~ x / 2 is approached. In- 
deed, the graph verifies that the theoretical exponent 
olb = 1/2 is adequately reproduced in the late time inter- 
val 10 3 < t < 10 4 . We measured the corresponding expo- 
nent a a for the mean-square A particle displacement in 
this regime as well, and found that the simulation results 
agree nicely with the mean- field exponent a a = 1/2, to 
the precision obtained in our simulations. 

We may also infer D sa 0.69 by matching our numer- 
ical integral of pB(t) with (x(t) A ), though incomplete 
knowledge of ps(i), particularly in the transient regime, 
introduces some errors to this estimate. A value D > 0.5 
indicates that more than one A particle is hopping to 
the same B particle site on average. When some A-B 
correlations develop there may be a higher density of A 
particles in the vicinity of a given B, which would in turn 
enhance the average A diffusion rate. 

Figure|21demonstrates that the predicted Gaussian dis- 
tribution Eq. l|15fl is not observed, despite the agreement 
in the time dependence of the mean-square displacement 
resulting from the distribution. Compared to a Gaussian 
with matching second moment, there is a distinct excess 
of essentially localized particles (with very small displace- 
ments), which necessarily implies longer 'tails' at large 
displacement values. A closer examination reveals that 
the simulation results for the A particle displacement dis- 
tribution agree over a large range of displacements with 
the normalized exponential distribution 



n A {x,t) 



\x\/L(t) 



(28) 



2L{t) 

where L{t) denotes a time-dependent characteristic 
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FIG. 3: The measured distribution of the A particle dis- 
placements corresponding to the system shown in Fig. |2] (B 
pair annihilation in one dimension, /i = 0.5, measured at 
t — 10000). Both Gaussian (dashed) and exponential fits 
(full line), uniquely determined by normalization and fixing 
the second moment with the simulation data, are depicted, 
on linear and logarithmic (inset) scales. The exponential fit 
works very well, although the value at and perhaps the tails 
of the observed distribution appear smaller than in the fit. 



length scale. As will be discussed below, the distribu- 
tion H28fl obeys dynamic scaling, so we determined L by 
matching the second moment of this function with the 
simulation data at t = 10000. As is evident from Fig. |3J 
the only significant deviations from this fit appear for 
particles with zero displacement and at the tails of the 
distribution, where simulation data contain greater error. 

We can understand the qualitative features of this non- 
Gaussian distribution as a result of the low connectivity 
of a one-dimensional lattice. In small regions where all 
B particles have been annihilated, the A species are tem- 
porarily localized and will have few subsequent chances 
to escape, thus lowering their effective diffusivity. Direct 
measurement of the B pair correlation function, as well 
as the slow B density decay clearly indicate that particle 
anti-correlations have developed. These anti-correlations 
enhance the probability of finding substantial regions of 
zero B particle density. However, as mentioned before 
in the discussion of the variation of the B particle ran- 
dom walk probability in Sec. II V Al a significant part of 
A movements probably result from hopping along with 
a particular B particle through several time steps (espe- 
cially in d = 1). Thus the effective diffusion coefficient 
is (at least temporarily) much higher for such particles, 
and this effect contributes results in the longer 'tails' in 
the displacement distribution. We may think of a local 
diffusion rate D{x,t) which is proportional to the distri- 
bution of available sites b{0, t) . Thus a particular particle 
will be subject to a temporally varying diffusion rate that 
depends on its trajectory through the field b(x,t), recall 



Eqs. I|12fl and l|13fl . just as the diffusing particle's tra- 
jectory must be considered for understanding diffusion 
on a static fractal. The distribution depicted in Fig. 
suggests a sort of phase separation into different pop- 
ulations: many particles are primarily subject to small 
diffusion rates, i.e., are mainly localized in regions with 
low B density, while some others experience quite large 
diffusivities (by remaining in areas of high local B parti- 
cle densities). 

The variation of local A diffusion rates is the result 
of inhomogeneities in b{x 1 t) and the coupling that al- 
lows an A particle to be carried along by a particular 
wandering B particle. Yet our earlier mean- field descrip- 
tion for the A species displacement distribution assumed 
that on average each A particle evolves with an identi- 
cal effective diffusion coefficient; i.e., the total number 
of hops for each A particle should be about the same. 
However, spatial inhomogeneities over scales larger than 
the distances traveled by typical A particles will cause 
the effective diffusivities for a particular time interval to 
vary spatially in a noticeable amount. This effect should 
be especially pronounced in the case of the B pair an- 
nihilation process because the emerging anti-correlations 
will leave many regions of space rather devoid of B parti- 
cles. A particles with larger effective diffusion coefficients 
probably follow a particular B particle for some time, 
since the B anti-correlations render them unlikely to be 
present in a region of high local B density. This distri- 
bution of B density spatial fluctuations should persist to 
a large extent throughout the simulation runs, which in 
turn should yield large variations in the average diffusion 
coefficient experienced by the A species. 

We may construct a simple model to accommodate 
such variations as follows. Suppose the number of B 
particles, i.e., the number of available hopping sites for 
the A species, is decreasing proportional to some nega- 
tive power of time, as is on average often the case in our 
simulations. If Eq. holds, we find for the fraction / of 
B particles that disappear during a small time interval 
At <C t: f(t) = olb At ft. Now assume that initially all of 
the A particles are diffusing ordinarily with an effective 
diffusivity D, and then after each small time interval At, 
a fraction / of the active A particles becomes localized 
and remains immobile for the remainder of the simula- 
tion. This suggests that the A particle distribution will 
be a superposition of normalized Gaussians of increasing 
width but multiplied by a factor to be found recursively 
from the condition that the A particle number be con- 
served. More precisely, in d = 1 this distribution becomes 
(with initial time to an d final time t = tu, M > 1): 

M 1 / 2 \ 

^■ t )-E^)(555P"»(-lk) (29) 

The prefactors p(tk) are then determined by keeping the 
integral of P(x,t) normalized to unity: ^2^—np(tk) = 1- 
Recursively, one thus arrives at p(to) = f(to), p{tk) = 
n^=o [! - for 1 < A: < M - 1, and at last 
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FIG. 4: The results from the 'continuous localization' model 
superimposed with the measured distribution of the A particle 
displacements corresponding to the system shown in Fig. [5] 
The model reproduces the faster than exponential drop-off 
in the tails of the distribution, but still underestimates the 
fraction of highly localized particles. The inset depicts the 
power law increase of {x{t) 2 A ) as computed from the model. 
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FIG. 5: The B particle density and corresponding A species 
mean-square displacement for the pair annihilation reaction 
2 B — > with reaction probability fi = 1 in d = 2. At long 
times, both plots indicate the expected logarithmic correc- 
tions to the mean-field scaling laws. 



p{tNi) = 1 [1 — f(tj)]. This picture assumes that the 
entire region surrounding an A particle is depleted of B 
particles nearly simultaneously, and that this region is 
not subsequently visited by other B particles. While this 
may be a decent approximation under the condition of B 
particle anti-correlations, it is rather more difficult to jus- 
tify in cases of emerging positive correlations (as we shall 
discuss below when offspring reactions are introduced). 
However, even then some B clusters are eventually elim- 
inated, leaving the nearby A particles localized, at least 
temporarily Furthermore the active A particles are dif- 
fusing normally with (x(t) 2 A ) ~ t in their local B cluster 
until they become trapped. 

Applying this simplified model to the pair annihilation 
system, we set as = 0.5, and we have chosen D = 0.2. 
We then summed Gaussian distributions for simulation 
times ranging from t — 5 to 100, localizing a fraction 
f(t) = as At/t at each integer t. The result is depicted 
in Fig. 01 after appropriately rescaling the axes to match 
normalization and the measured second moment of the 
simulation data. Apart from large deviations for small 
displacements, this 'continuous localization' fit appears 
to agree remarkably well with the simulation results. In- 
deed this generated A displacement distribution even be- 
gins to decrease more quickly near its tails, as may also 
be observed in the simulation data. The disagreement 
between the data and model may of course be traced to 
the simplifications assumed in its construction. However, 
the error also might arise from choosing the initial dis- 
tribution: at the beginning of the simulation the number 
of B particles decreases sharply before anti-correlations 
have developed, and the above estimate for the fraction 
of localized A particles is not valid for At/t ~ 1. Despite 



its shortcomings, the model seems to capture most of the 
features we have observed in the A species displacement 
distribution. We also computed the time dependence of 
the mean-square displacement using the 'continuous lo- 
calization' model. Following a transient behavior, the in- 
set in Fig.^Jshows (x(t) A ) ~ t 55 , close to the expected 
power law with 1 — a a = 0.5 (which we also observed in 
the simulation) . By trying a variety of as values one can 
generate different approximate power laws, but the mea- 
sured subdiffusive A displacement exponent was always 
found to be slightly larger than 1 — olb- 

To examine the non-Gaussian character of tia{%, t) 
further, we measured its higher moments as a function 
of time. In d = 1, Eq. (|17fl yields for the even mo- 
ments of the mean-field Gaussian distribution: (x(t) A k ) = 
{2k — 1)!! (x(t) A ) k . While this factorization of higher 
moments still holds to the accuracy of our simulation 
data, we found the prefactors = (x(t) A k )/(x(t) A ) k 
in our measurements to differ from the predicted val- 
ues C2 = 3, C3 = 15, and C4 = 105: we measured 
the considerably larger values C2 ~ 5.5, C3 « 55, and 
C4 k 1100. These numbers are actually closer to the 
values one would obtain from the exponential distribu- 
tion (gSJl, namely c k = {2k)\/2 k , i.e., c 2 = 6, c 3 = 90, 
and C4 = 2520, but still off by a factor of about two 
for the higher moments. Thus the exponential fit cannot 
entirely describe the measured distribution either. Nev- 
ertheless, the factorization property (x(t) 2 A ) oc (x(t) 2 A ) k 
found in the simulation data is significant, for it indi- 
cates dynamic scaling: once the time dependence of the 
characteristic length scale L(t) = (x(t) 2 A ) 1/2 ~ t {1 - aA)/2 
is factored out, the shape of the probability distribution 
should remain constant in time. The distributions l|15|) 
and (|28|l clearly display this feature. Since the proba- 
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FIG. 6: The measured distribution of the final A particle 
displacements corresponding to the system shown in Fig. |S] 
(B pair annihilation in two dimensions, fj, = 1), shown with 
Gaussian and unnormalized exponential fits. 



bility distribution is fully determined by all its moments, 
our measurements of the first four moments indicate that 
the true ua(x, t) obeys dynamic scaling as well. 

At g? c (2) = 2 one obtains the typical logarithmic cor- 
rections to the mean- field result, viz. Eq. © for the B 
particle density, and Eq. i|21|) for the A species mean- 
square displacement. Simulations were also carried out 
at various values of /x for a two-dimensional system. Fig- 
ure [S] illustrates the apparent agreement with the above 
predictions at long times for simulations run with anni- 
hilation probability /i = 1 (averaged over 5 runs on a 
100 x 100 lattice). However, such a large annihilation 
rate severely suppresses the A particle diffusion, so one 
should not really draw too firm conclusions from these 
measurements. We have again measured the correspond- 
ing distribution of the A particle displacements. In Fig. [HI 
we see that nA(x,t) is non-Gaussian and resembles its 
one-dimensional counterpart (Fig. though the devia- 
tions from Eq. (|15fl are less pronounced. When // <C 1, 
anti-correlations develop over much larger time scales and 
as the corrections are only logarithmic, for fi = 0.01 (20 
runs on a 200 x 200 lattice) we merely recovered the mean- 
field results in the regime accessible to our simulations. 
Though not shown here, the B density decay was cap- 
tured by the mean-field power law on the time scales of 
our runs, and the A particle mean-square displacement 
followed the integral of that quantity after some initial 
transient behavior. Accordingly, the A particle distribu- 
tion was essentially Gaussian in this situation. 



C. Triplet annihilation 35^0 

Next we consider the triplet annihilation reaction 
3 B 9. The upper critical dimension here is d c (3) = 1. 
Hence in one dimension we expect logarithmic corrections 
to the mean-field power law B density decay, as given in 
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FIG. 7: The B particle density and corresponding A species 
mean-square displacement for the triplet annihilation reaction 
3B — » with j« = 1 in d = 1, both displaying the expected 
logarithmic corrections to the mean-field scaling laws. 

Eq. 10. According to the simple mean-field picture, this 
leads to Eq. (|22|l for the A species mean-square displace- 
ment. Figure [7] shows our simulation results (averaged 
over 20 runs on a lattice with 10000 sites, annihilation 
probability fj, — 1). We are able to clearly detect the log- 
arithmic corrections even at fj, = 1 since this reaction is 
a much slower process than pair annihilation (requiring 
three particles to meet on a site) . The power law regres- 
sion on {x(t) 2 A ) I '(hit) 1 / 2 yields a value of 1 — oca ~ 0.52, 
whereas the expected value is 0.50. Supposing the mean- 
field result fully applies here, this gives an idea of the 
overall precision of our simulation data. The A particle 
displacement distribution also agrees well with the pre- 
dicted Gaussian, apart from a slight excess of presumably 
localized particles around (x A ) = and correspondingly 
longer 'tails' of the distribution. Yet the deviation is 
much smaller than in the pair annihilation case because 
the anti-correlations are less pronounced here. 

For d = 2 > d c , the mean-field result j£| should pro- 
vide a correct description, i.e., for n = 3 in the long-time 
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FIG. 8: The B particle density and A species mean-square 
displacement for the triplet annihilation reaction 3 B — > 
(fj, = 1) in d = 2, confirming the mean-field power laws. 
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limit t » r 3 : p B {t) ~ i" 1 / 2 and (x(t) 2 A ) ~ t 1 / 2 . Fig- 
ure demonstrates agreement with these mean-field pre- 
dictions after 10 3 or fewer timesteps (from 3 runs on a 
100 x 100 lattice with again p = 1). In fact, to our res- 
olution the mean-square displacement of the A species 
converges to the mean-field power law markedly faster 
than the B particle density. As yet we have no explana- 
tion for this surprising observation. 

D. Quartic annihilation 4B — > 

The kinetics of quartic annihilation should be aptly de- 
scribed by the mean-field rate equation Q in all physical 
dimensions since d c (4) = 2/3. Setting n = 4 in Eq. (JSJ) we 
thus expect p B (t) ~ t^ 1 ^ for sufficiently large t t&. 
According to Eqs. (E) and 0, (x(t) 2 A ) ~ t 2 / 3 . We 
ran simulations in two dimensions on a 100 x 100 lattice 
(with (i=l, averaged over 8 runs) and found good agree- 
ment with these mean-field scaling laws at sufficiently 
long times both for the B particle density decay and 
the A species mean-square displacement. For the anni- 
hilation processes, larger n values imply longer crossover 
time scales r„ since n particles must meet for a reaction 
to occur, see Eq. Indeed, by comparing with the re- 
sults for the triplet reactions, we noticed that the time 
interval of transient behavior prior to convergence to the 
asymptotic mean-field power law was at least an order of 
magnitude longer for the n = 4 system. As in the triplet 
simulations, the convergence to mean-field behavior oc- 
cured faster for the A particle mean-square displacement 
than for the total B density. 



V. ACTIVE TO ABSORBING STATE PHASE 
TRANSITION AND LOCALIZATION 

After this preliminary study with pure B annihilation 
kinetics, we now turn our attention to the primary goal 
of our investigation, namely anomalous diffusion on a dy- 
namic fractal (see Fig. . Each reaction discussed below 
exhibits a continuous phase transition separating active 
and absorbing stationary states. At the critical point, the 
B particle distribution is known to be fractal in space and 
time. Thus the assumption of a homogeneous B particle 
distribution is significantly violated, and we expect the 
mean-field description for the A species propagation to 
be inadequate. Our aim has been to measure and char- 
acterize the deviations from the mean-field predictions. 
All of the phase transitions in the B particle system to 
be discussed here are described either by the directed 
percolation (DP) or parity conserving (PC) universality 
classes. The accepted critical DP exponents are listed in 
Tab.HJ As mentioned in the Sec. HI a PC phase transition 
is observed in d — 1 when A = in the reaction scheme 
(JTJ , and m is even; in higher dimensions the (degenerate) 
critical point occurs at zero branching rate and is gov- 
erned by mean-field exponents 0. The PC exponents in 



d = 1 are also given in Tab. [I] 

Knowledge of the B particle density behavior of these 
systems in the active and absorbing states near the phase 
transition is also important, both in locating the critical 
point and for comparison with the asymptotic critical 
scaling laws. In fact, in simulations we will inevitably 
always be slightly away from the precise critical control 
parameter values. However, in the vicinity of the phase 
transition we expect the system to exhibit the critical 
power laws for some time interval before crossing over 
to sub- or supercritical behavior. The closer the system 
parameter values are to the critical ones, the longer this 
critical regime lasts, as also suggested by Eq. (|27J) . and 
the more precisely critical exponents can be measured. 
When m is odd (independent of A), we expect systems 
in the absorbing state (i.e., dominated by annihilation 
reactions) to exhibit behavior dictated by Eq. ©, the 
solution to B particle radioactive decay, with some ef- 
fective rate A e ff, dependent on A, a, and p. To be more 
precise, let us consider the reactions Q with m — 1 and 
n = 2 as is the case in many of the situations examined 
below. Recalling that the branching reactions are local 
(offspring particles are placed on the parents' neighboring 
sites), in low dimensions there is a significant probability 
that the parent and offspring will meet in the next few 
time steps after the branching reaction and undergo pair 
annihilation. Together, the branching and subsequent 
annihilation reactions generate the decay B — > 0, even if 
A = on the outset. Generally, this is true when n = 2 
and m is odd. However, when n = 2 and m is even, local 
parity conservation eliminates the possibility to generate 
'spontaneous' death processes. Thus in the PC universal- 
ity class the inactive phase is characterized by the power 
laws of the pair annihilation reactions 2B — > 0, whence 
in subcritical systems one should asymptotically observe 
the behavior given in Eq. Q with d = 1. In the ac- 
tive phases of both DP and PC systems, the B density 
will approach its stationary value exponentially. This 
follows immediately from linearizing the corresponding 
mean-field rate equations. 



A. Directed percolation universality class, d = 1 

To search for the DP transition in one dimension, we 
considered the reactions with m = 1 and n = 2, with 
A = 0.01 and p = 1 fixed while a was varied. Notice 
that setting p — 1 eliminates the possibility of multiple 
B particles occupying the same site, so that the B par- 
ticle density exactly corresponds to the density of avail- 
able lattice sites for A particle hopping. According to 
the list in Tab. H p B (t) ~ t~ 0159 at criticality, which 
implies within the mean- field approximation (a a = ot B ) 
that {x(t) 2 A ) ~ t - 841 . Our closest estimate of the critical 
point is er c 0.8975. Figure [5] shows the power law de- 
pendence of p B {t) and (x(t) 2 A ), obtained from 26 runs on 
a lattice with 10000 sites. The double-logarithmic linear 
regressions yield values a B — 0.161 and <xa — 0.167 over 
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FIG. 9: The B particle density and A species mean-square 
displacement for the DP reactions B^0(A = O.O1),B^2S 
(a = 0.8975), 2B -> (p = 1) in d = 1. 



more than three decades, i.e., a a ~ «s ~ 0.006. As we 
estimate the numerical uncertainty of the measured ex- 
ponents to be at least ±0.01, we observe agreement both 
between cub and the DP prediction as well as between as 
and a a within our error bars. We also computed local ex- 
ponents, as defined in Eq. (|2*oT) setting b = 2. Figure ITUI 
indicates that after following some transient behavior, 
oiA ~ 0.158, in excellent agreement with the critical DP 
value. Thus we conclude that our mean-field prediction 
of the A particle mean-square displacement time depen- 
dence works excellently for this reaction at least on the 
time scales we were able to access with our simulations. 

In Fig. ^2 (inset) we demonstrate the algebraic de- 
pendence of the B particle correlation function on par- 
ticle separation, indicating the spatially fractal structure 
at criticality. The data were taken at t = 50000, aver- 
aged over 9 simulations. The measured exponent 2j3/v 
agrees well with the DP prediction 0.50 in d = 1 over two 
decades. We remark that for small distances we observed 
that the correlation function scales as Cs(|x|) ~ |a;| _/3 / I/ , 
i.e., with precisely one half the asymptotic scaling ex- 
ponent. The origin of this can be understood as fol- 
lows: at low B particle densities, the site occupation 
number n(x) is either zero or one; hence for \x\ < £ 
(n(x)n(0)) ~ (n 2 ) = (n) ~ ps ~ S y ~^^ v . The propor- 
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FIG. 10: The local exponent oa corresponding to data shown 
in Fig. (with 6 = 2) and Fig. EU (with b = 4). 




FIG. 11: The A and B particle pair correlation functions 
Ca(x) and Cb(x) (inset) at the DP critical point (d = 1, 
measured at t = 50000). Cb(x) ~ \x\~ 213 ^ is approximately 
scale-invariant. The A particles are very likely to be found on 
the same site, but there is a slight negative correlation with 
nearby sites. Ca(x) is essentially zero at distances \x\ > 10. 



tionality factor here must be a scaling function /(r/£), 
and demanding that the £ dependence must cancel at 
criticality, we obtain the aforementioned scaling law. 
As we are not precisely at the critical point, we see a 
crossover to exponential decay at large distances. 

FigureEldepicts the measured full A particle displace- 
ment distribution at t = 50000. It displays an excess of 
localized particles and corresponding long tails of the dis- 
tribution, similar to the pair annihilation case. But the 
deviation from a Gaussian distribution appears much less 
pronounced in this system. The power law dependence 
of B particle spatial correlations qualitatively suggests 
the observed distribution. Such strong correlations yield 
macroscopic regions of relatively large B densities which 
an A particle may traverse and thus acquire a large dis- 
placement. However, such clusters necessarily indicate 
compensating large regions of very low B particle den- 
sity, where the A particles become highly localized. The 
A particles at the fractal 'boundary', which by its nature 
accounts for a significant fraction of the system volume, 
appear to behave in accord with the mean-field predic- 
tion. This is suggested by the fact that the data in Fig. 1121 
agree well with a Gaussian distribution for intermediate 
displacement values. 

Yet it was in these fractal boundaries that we expected 
the deviation from the mean-field approximation to arise. 
One possibility we had imagined is that A particles in 
the boundary region are driven towards denser B regions 
owing to the particle density gradient expected from the 
formation of clusters, thus destroying the spatial homo- 
geneity of the A particles and resulting in a higher net 
diffusion rate than assumed by mean-field theory. This 
behavior can be distinguished by computing the A- A par- 
ticle correlations. If the described large-scale A-B cou- 
pling were strong enough, then we would expect to see 
some A particle correlations as they congregate in regions 
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FIG. 12: The measured distribution of the A particle dis- 
placements corresponding to the system shown in Fig. [5] 
(DP critical point for the B system in d = 1, measured at 
t = 50000). 



of high B density. But these correlations were not ob- 
served, as demonstrated in Fig. ^2 We may attribute 
the inability of the B particles to induce significant cor- 
relation in the A system to the low connectivity of a 
one-dimensional lattice. However, Ca(0) does indicate a 
significantly increased probability to find multiple occu- 
pation of a given site. A weak negative correlation seems 
to have developed for < \x\ < 5. Such a distribution 
may be the result of the 'piggy-back' effect discussed ear- 
lier: several A particles may in fact follow a single B par- 
ticle through a number of time steps. If the B particle 
is subsequently annihilated, the 'piggy-backing' A parti- 
cles are all (temporarily) localized at their current site. 
The negative correlation may just be the compensation 
for the effective mutual attraction of a group of nearby 
A particles all following the same B. 

We have also measured the higher moments of the A 
displacement distribution nA{x,t). We found (x{t)\) = 
Cfc (x(t) 2 A ) k , indicating dynamic scaling, albeit with C2 = 
3.3, C3 = 20, and C4 = 185 compared with the values 3, 
15, and 105 that would result from a Gaussian. A few 
distinctions between the pair annihilation case and this 
DP phase transition may account for the significant dif- 
ferences in the corresponding deviations from the mean- 
field A displacement distribution. Firstly, the B particle 
density decays much slower in the DP case (as w 0.16 as 
compared with 0.5). This slower decay may cause the DP 
distribution to be dominated by the 'active' A particles. 
For a fixed system size we will also observe the DP system 
for longer durations. Since the B particles undergo ordi- 
nary diffusion while reacting, these extended time scales 
imply that previously localized regions (where the local 
B density is zero) are likely to be visited by wandering 
B particles. Thus the localized portion of the A particle 
distribution becomes smeared as these A particles are re- 
activated, which makes it more likely that all A particles 
experience roughly the same average effective diffusion 
coefficient throughout the simulation (as assumed in the 



mean- field approach). Finally, at the DP phase transi- 
tion, the B particles are positively correlated, whereas 
the pair annihilation process induces negative correla- 
tions. Positive correlations may also contrive to weaken 
particle localization, since the A species in regions of high 
local B density are less dependent on single B particles 
for their mobility. We also note that as a consequence of 
the effective 'slaving' of the A species by the B particles, 
the correlation length exponents v and v t = zv should 
be identical for both the active to absorbing transition in 
the B system and the induced A localization transition. 
Indeed, to the accuracy we could determine those expo- 
nents by means of Eq. I|27l) and the method described in 
Ref. [T3, this equality appeared to hold. 

The expected time independence of the equal-time B 
density correlation function and the extent of the B par- 
ticle clusters explain the time independence of the A dis- 
placement distribution. Recall that the origin of the non- 
Gaussian distribution is that the A particles experience 
different effective diffusivities depending on their loca- 
tion with respect to B clusters. What determines an 
A particle's placement in the distribution is the average 
diffusion coefficient experienced over the length of the 
simulation, neglecting the possibility of biased diffusion 
due to B particle gradients. For the shape of the dis- 
placement distribution to remain constant in time, the 
shape of the distribution of time-averaged effective dif- 
fusion rates most likely remains constant as well (even 
while on average the effective A species diffusion rate 
is decreasing with time). Subsequent measurement of 
the A diffusivity distribution supports this interpreta- 
tion, see Fig. [57| in Sec. IVII below. This is only possible 
if A particles are able to remain in a dense B region 
for a large portion of the simulation, made possible by 
large cluster sizes implied by the power law correlation 
function, while others remain in a low density region for 
long durations. These two persisting extremes maintain 
the non-Gaussian nature of the A species displacement 
distribution. However, A particles at the fractal bound- 
ary, which comprises a considerable volume of the system 
and thus contains a large fraction of the A particles, may 
all see roughly the same effective diffusion rate over the 
length of the simulation, which results in the Gaussian 
mid-region of nA(x,t). Furthermore, the stationarity of 
the shape of the distribution of effective A diffusivities is 
likely the reason we see such good agreement between a a 
and as- For instance, even if the global average instan- 
taneous diffusion rate evolved as D(pB(t)) (as we should 
expect), a changing distribution shape, such as increasing 
enhancement of the phase separation between 'localized' 
and 'active' regions while maintaining the proper global 
diffusion rate, would induce additional time dependence 
and cause deviations from Eq. I|18fl . However, we have 
no precise argument for the apparent time-independence 
of these distributions. 

We also sought to verify the predicted behavior away 
from the phase transition. Setting a — 0.910, sufficiently 
above the critical branching rate, we observed almost im- 
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FIG. 13: A species mean-square displacement in the super- 
and subcritical regimes for the DP processes of Fig.|5|(d = 1). 
For a = 0.910 > a c , one finds normal A diffusion. In the 
subcritical regime (er = 0.890), the B particle density de- 
cays exponentially with A c ff » 0.0002, and the mean-square A 
displacement approaches the asymptotic value exponentially 
with the same time constant. (Both subcritical data sets are 
rescaled to a unit prefactor.) 



mediate convergence of the B particle density to its ac- 
tive state saturation value p s — (pB(t)) ~ 0.26, with 
standard deviation 0.009. As expected and depicted in 
the inset of Fig. ^3 the A species then exhibits normal 
diffusion (a a — 0). (The slight deviation at the end of 
the single run on 10000 sites may be a finite-size effect.) 
By Eqs. I|18l) and l|16l) we also expect the prefactor to 
be 2Dp s , wherefrom we estimate D rj 0.39. For ordinary 
diffusion, in the absence of any correlations, D — 0.5. In- 
deed, our algorithm dictates choosing a hopping direction 
first, and then checking for the presence of a B particle at 
the destination site. Thus on average one of two neigh- 
boring A particles will hop to a B site. To observe sub- 
critical behavior, we set a = 0.890. Figure H"3l illustrates 
the convergence to the expected exponential density de- 
cay with an effective decay rate A e g ~ 0.0002 (10 runs 
on a system with 10000 sites). Initially (for t < 1000), 
a brief quasi-critical regime is visible with power-law de- 
cay ps(t) ~ t~ - 26 and correspondingly (x(t) A ) ~ t ' 77 . 
For t > 2000 we then observe exponential convergence to 
Pb = and the stationary value for (x(t) A ), both with 
the same time constant A e ff w 0.0002. 

As discussed above, we expect the macroscopic proper- 
ties of the DP phase transition to be independent of the 
value of A, for d < 2 even persisting to A = 0, i.e., the ab- 
sence of spontaneous decay. To check this, we ran simula- 
tions for the reactions with m = 1 and n — 2, setting 
A = 0, fj, = 1, and varying a. We estimated a c « 0.8930, 
slightly below the value found when A = 0.01. Figure ITU 
demonstrates the critical behavior deduced from data 
taken over more than three decades (20 runs in a lat- 
tice with 20000 sites), before the transition to subcritical 
behavior becomes evident. We found as ~ 0.158, in su- 
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FIG. 14: The B particle density and A species mean-square 
displacement for the DP reactions (branching and annihi- 
lating random walks with odd offspring number) B — > 2 B 
{a = 0.8930) and IB -> (jj, = 1) in d = 1. For t < 200000, 
the critical power laws are observed. Subsequently a crossover 
to subcritical behavior can be seen. 

perb agreement with the expected DP value 0.159. We 
measured a. a — 0.157, so as — a a ~ 0.001. Again any 
significant deviations from the mean-field prediction at 
the critical point, if present at all, must arise at time 
scales inaccessible to our simulations. 

However, the local exponents a a and as, setting 6 = 4 
in Eq. I|2tj|) . turned out not to be constant. In the regime 
from which we inferred the global values of as, the local 
values oscillate about their averages by ss 0.01, whereas 
the local a.A is a much smoother function in time (see 
Fig. ITUI). The steady increase in a a (and as) indicates 
that the system is actually in the inactive phase. We 
found the increase in the local as exponent to be signif- 
icantly faster. However, this is not an indication of devi- 
ation from the mean-field prediction for the mean-square 
displacement of the A species. We in fact computed nu- 
merically the integral of the B particle density and com- 
pared it with the measured mean-square A displacement. 
We found that for t > 10000, the error between the two 
measurements was less than 1%. We estimated D — 0.39 
for this reaction in order to obtain the best match be- 
tween the measured mean-square A displacement and the 
integral of /Os(i). This value is in good agreement with 
estimates from other reactions. Both the B and A parti- 
cle correlation functions and the A species displacement 
distribution were basically identical to Figs. ITT1 and Ir21 



B. Directed percolation universality class, d = 2 

To investigate the properties at a DP transition in two 
dimensions, we used the same reactions and rates as in 
d = 1, namely B -> (A = 0.01), B -> 2B (varying 
cr), and 2B -> (p = 1). For d = 2, Tab. [Q predicts 
PB(t) ~ t~ 0A51 , which implies, according to our mean- 
field description, that (x(t) 2 A ) ~ t - 549 . Our best estimate 
of the critical branching rate is er c « 0.2233. Figure IT51 
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FIG. 15: The B particle density and A species mean-square 
displacement for the DP reactions B^0(A = O.O1),B^2B 
(cr = 0.2233), 2 B -> (jj, = 1) in d = 2. 



shows the power law dependence of ps(t) and (x(t) A ), 
inferred from 40 runs on a 800 x 800 square lattice. The 
initial A density was set to 0.01 (as opposed to the usual 
0.5) here. The power law fits imply D — 0.256, very 
close to the expected 0.25 from ordinary diffusion in two 
dimensions (since on average one A particle will hop to 
each available B site) . The slight error in computing D is 
probably a result of the transient behavior at the begin- 
ning of the simulation, where the B density is still larger 
than the power law fit would predict. Figure [Tol shows 
the local exponent a a (computed with b — 2), indicating 
that the system is indeed subcritical, as a.A is increasing 
with time for large t. But the minimum plateau value is 
very close to the predicted qa 0.45. Figure El demon- 
strates that the B particle distribution is still fractal at 
t = 50000, yielding an effective exponent value 2(3 jv close 
to 1.59 (as expected from DP) at intermediate distances. 
Figure ED depicts the A species displacement distribu- 
tion at t = 50000 together with the mean-field Gaussian 
and an (unnormalized) exponential fit of the distribution 
tails. As in previous cases, the 'tails' of the distribution 
are longer than a Gaussian would suggest, and instead 
appear to be matched well by an exponential. 



0.5 -i 
0.49 - 
0.48 
0.47 
0.46 
0.45 




10000 f 15000 20000 25000 



FIG. 17: The B particle pair correlation function Cb(%) at 
the DP critical point (d — 2, measured at t = 50000). 



At last, we verify the DP transition in two dimensions 
when A = 0, and compare the A particle anomalous 
diffusion to the nonzero A case, using the same reac- 
tion rates as before in d = 1. Because as is so large 
for DP in d = 2, we had to use quite sizeable systems 
(800 x 800) to measure exponents over several decades in 
search of the critical point. Thus our determination of 
cr c was not as precise. For the B particle density and A 
species mean-square displacement as functions of time, 
the measured exponents are as ~ 0.47, approximately 
0.02 larger than the DP value, and cva ~ 0.48. We found 
our system close to but below criticality, and therefore 
observed time-dependent local exponents a a and ag. We 
also completed fewer (20) runs for this system, and so 
expect a larger error on our measurements of the expo- 
nents. The power law fit prefactors imply D w 0.29, in 
fair agreement with our expectation of 0.25 and previ- 
ous measurements. In summary, we have not uncovered 
any significant deviations from the behavior observed for 




FIG. 16: The local exponent oa corresponding to data shown 
in Fig. [Tol ( with 6 = 2). 



FIG. 18: The measured distribution of the A particle dis- 
placements corresponding to the system shown in Fig. 1151 
(DP critical point for the B system in d = 2, measured at 
t = 50000). The dashed line indicates a Gaussian fit, and the 
solid line is represents an unnormalized exponential fit to the 
distribution tails. 
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FIG. 19: The B particle density and A species mean-square 

displacement for the PC reactions B -> 3B (a = 0.2175), FIG. 21: The B particle density and A species mean-square 
2B — > (fj, = 0.5) in d = 1. displacement for the PC reactions B — > 5 B (a = 0.2795), 

2B -> (/i = 1) in d = 1. 



A > 0. 



C. Parity conserving universality class, d = 1 

In order to observe the phase transitions from active to 
inactive/ absorbing states in the PC universality class, we 
looked at branching and annihilating random walks with 
even number of offspring particles, i.e., set A = 0, n — 2, 
and either m = 2 or 4 in the reaction scheme Q. In the 
former case, 2B — > combined with B — > 3B, we were 
forced to set the annihilation probability to /i = 0.5 in or- 
der to detect the phase transition (at varying branching 
probability a). FigurelT9lshows the power laws for the B 
particle density decay and the A species mean-square dis- 
placement as a function of time near the phase transition 
(at cr c « 0.2175). Our measurement of as = 0.269 agrees 
fairly well with the expected PC exponent 0.285 (from 20 
runs, each with 60000 sites). Furthermore, a a = 0.271 
is in excellent agreement with as- The power law fits 
shown imply D — 0.45 for this reaction. Figure [5D| de- 
picts a measurement of the B density correlation function 
at t = 500000, and the predicted PC exponent ratio 1.15 
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FIG. 20: The B particle pair correlation function Cb(x) at 
the PC critical point in d = 1 (measured at t — 500000). 



according to Tab.|U Measurements of Cb (x) at t = 50000 
and t — 100000 yielded distributions very similar to that 
shown in Fig. [201 Figure \\\in Sec. [i] shows a space-time 
plot for these reactions illustrating the fractal nature of 
the process at criticality. 

We also observed the critical behavior for the PC reac- 
tions 2B — » (/i = 1) combined with B — > 5B. Setting 
/i = 1 helps to minimize the number of necessary random 
numbers, as well as prohibits the multiple occupation of 
lattice sites and thus avoids discrepancies between the B 
particle density and the density of available sites for the A 
species. In Fig.|2]we display the B density and A species 
mean-square displacement as a function of time near the 
phase transition (at a c ~ 0.2795), implying D = 0.50 
(averaged over 20 runs on the lattice with 60000 sites). 
Indeed we again find good agreement as ~ 0.274 « a a, 
as well as with the simulation data depicted in Fig. 1191 
and the predicted PC value. In addition, we measured 
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FIG. 22: The measured distribution of the A particle dis- 
placements corresponding to the system shown in Fig. l21l (PC 
critical point for the B system, d = 1) at different times as 
indicated. 
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FIG. 23: The paths of the five A particles with largest dis- 
placements in a simulation contributing to the data plotted 
in Fig. 1211 Notice the long intervals of particle localization 
interspersed with brief periods of high mobility. 
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the A displacement distribution at three distinct times: 
t = 50000, t = 100000, and t = 600000. By scaling 
the ordinate by a factor s — (ta/ti)^ -0 " 4 " 2 , and the ab- 
scissa by 1/s, the three curves are seen to be essentially 
identical in shape, see Fig. |221 which supports the dy- 
namic scaling conjecture that the A particle distribution 
maintains its shape as it evolves in time. In Fig. 1231 we 
show the paths of five A particles in one of the simula- 
tions contributing to the data in Fig. |^ The selected A 
particles are those with the greatest displacements. The 
trajectories consist of periods of localization (very little 
activity) interspersed with periods of much movement. 
This diagram resembles those for other cases of anoma- 
lous diffusion such as Levy flights. We also examined a 
subcritical system, setting a = 0.277. Figure [32J shows 
an initial critical regime with pB(t) ~ t , crossing 
over very quickly to the subcritical behavior dominated 
by pair annihilation, ~ t~~ ' 5 . We see that through- 

out these regimes the mean-square displacement of the A 
particles agrees remarkably well with the integral of the 
B density, setting D = 0.425 (from 7 runs, each with 
40000 sites). 

To uncover the origins of the non-Gaussian distribu- 
tions observed in the systems thus far examined, we also 
considered a simple model in which the B particle distri- 
bution remained uncorrelated, while still exhibiting the 
desired overall density decay. The algorithm allowed the 
B species to diffuse through the lattice without site exclu- 
sion, but removed sufficiently many of them at random 
so that the B density followed the prescribed global be- 
havior. For example, we forced the B density to decay 
as observed in Fig. 1211 near the PC phase transition in 
d = 1. As would be expected, the time dependence of the 
mean-square displacement of the A species agreed with 
the integral of the B particle density with the absolute 
error bounded by 10. The required choice for this agree- 
ment was D = 0.95, a puzzlingly large value, considering 
all previous observations gave D w 0.5. Since in this algo- 
rithm the B particles are not being continuously created 
and annihilated, the A-B coupling is actually stronger, 



FIG. 24: The B particle density and A species mean-square 
displacement for the PC system of Fig. 1211 in the subcritical 
regime (a = 0.277 < a c ). The crossover from the critical 
power laws to those of the inactive phase is clearly visible. 



and thus on average there is probably a higher density 
of A particles surrounding a typical B site, thereby in- 
creasing the effective A diffusivity. We may now compare 
the A displacement distributions at t = 100000 for the 
homogeneous and true PC systems, see Fig- ESI Different 
effective diffusivities required an area-preserving reseat- 
ing. While there does appear to be some deviation from a 
Gaussian even in the artificial model, the figures suggest 
that the deviation is not as great as in the true PC sys- 
tem with B particle correlations. Though not shown, a 
rescaled version of the radioactive decay A displacement 
distribution agrees well with the 'artificial' distribution 
shown here, as one might expect. We suggest that the 
low connectivity of the one-dimensional lattice necessar- 
ily causes significant deviations from the average effective 
diffusivity D(pg(t)) for the A species, which then yields 
a non-Gaussian distribution of the A particle displace- 
ments. However, Fig. |23 demonstrates that B particle 
correlations are also significant in producing deviations 
from the mean-field prediction. 

We also explored a few other algorithms to verify that 
our results were not particular to our implementation. 
First we adapted the A species hopping algorithm, so 
that whenever an A particle could hop, it would with 
probability 1. Previously, A particle hopping was imple- 
mented by choosing a direction to hop, and then checking 
to see if the destination nearest neighboring site was oc- 
cupied by a B particle. We executed the maximal A 
particle hopping algorithm for a = 0.2795, our best esti- 
mate of the critical point for the PC system with m = 4. 
The time dependence of the mean particle displacement 
agreed very closely with the integral of the B particle 
density, with D = 0.87. At t = 100000 the A displace- 
ment distribution coincided with that in Fig. 1221 apart 
from a scaling factor to account for the discrepancy in D. 
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FIG. 25: The measured distribution of the A particle dis- 
placements corresponding to the system shown in Fig, l21l (PC 
critical point for the B system, d — 1) , compared with the one 
resulting from the artificial model with uncorrelated B parti- 
cles. Notice that different effective diffusivities were applied 
to obtain equal second moments. 

The increase in this value simply indicates that A parti- 
cles will hop whenever a neighboring site is available. 

Next we evolved the B particles until t = 500, when 
the B density was down to 0.07. Thus clusters and rel- 
atively vacant regions should have begun to form. The 
algorithm then placed A particles on the existing B parti- 
cles and studied the A diffusion thereafter. Though after 
a longer transient, the mean-square A displacement ap- 
proached the integral of the B density (setting D = 0.5). 
Finally we examined an algorithm that removed A par- 
ticles from the system after they had been localized for a 
certain time interval (here chosen as 10000 Monte Carlo 
steps). We found the number of active A particles to 
be a sharply decreasing function of time. The localized 
particles being removed from the computation had lower 
effective diffusion rates than the remaining active parti- 
cles, overall increasing the mean-square A displacement. 
This observation supports the hypothesis that the effec- 
tive diffusion rate for the A species is not uniform, indi- 
cating that in the previous models, a significant fraction 
of the A species was localized over time scales at least as 
large as 10000 Monte Carlo steps. Thus most A particles 
are not truly diffusing through the dynamic fractal, but 
rather execute an occasional hop when passed over by B 
particles. 



VI. CONCLUDING REMARKS 

We have numerically investigated the reaction- 
controlled diffusion model introduced in Rcf. for vari- 
ous B species reaction systems, and studied the ensuing 
anomalous diffusion for the A particles on the emerging 
dynamic fractal B clusters. We found excellent agree- 
ment for all examined systems of the A species mean- 
square displacement (x(t) A ) with the integral of the mean 
B particle density (ps(i)}, at least within the accuracy of 
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FIG. 26: The measured A species displacement distributions 
for the majority of systems investigated here (PA and TA rep- 
resent the pure B particle pair and triplet annihilation pro- 
cesses, respectively). We see that most of the data collapse 
to roughly the same scaling function (except for the result 
for one- dimensional B pair annihilation with the strongest 
anti-correlations). The observed distribution clearly deviates 
from the mean-held Gaussian both at small and large dis- 
placements. 



our data. However, the mean-field Gaussian distribution 
(|15|l for the A particle displacements was not observed. 
We have argued that the primary factor in creating this 
deviation is a variation in the effective diffusion rate of 
the A particles, which is proportional to the number of 
hops executed by an A throughout the simulation run. 
The shape of the resulting A displacement distribution 
was seen to be the same at different times and also over 
many different systems (see Fig. 1261) . thus suggesting that 
the distribution shape of effective diffusion rates also re- 
mains the same over different time scales and systems. 

The contributing factors to producing the diffusion 
rate distribution were identified as persisting spatial fluc- 
tuations in the local B particle density (enhanced in sit- 
uations with strong B correlations), the low connectiv- 
ity of the one- and two-dimensional lattices examined 
here, and the fact that the A species tend to be carried 
along with diffusing B particles (the 'piggy-back' effect) . 
As mentioned above, the distribution of effective diffu- 
sivities is constrained in that the average diffusion rate 
should equal D(ps(t)). The evolution of this average 
diffusion rate according to t^ 1 J Q pB(t')dt' (see Fig. I27J) 
coupled with the time independence of the A distribu- 
tion shape, i.e., dynamic scaling, essentially yields the 
mean-field time behavior of (x^(t)), Eq. (|18|) . Also, as 
the mean-square displacement is an integral quantity, the 
noise associated with fluctuations in the B particle den- 
sity becomes suppressed. Consequently, under the as- 
sumption that a a — cub to a very good approximation 
at least, one may measure critical exponents of the re- 
acting B particle system via the passive A species using 
fewer simulation runs and with reduced statistical noise. 

Fluctuation effects really become manifest in the 
higher moments only, and in the overall shape of the A 
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FIG. 27: The distribution of effective A diffusivities, obtained 
as the average number of hops per unit time, at various times 
during simulations for critical DP processes for the B system 
in d = 1. The inset shows that the average diffusion coefficient 
becomes identical with the temporal average of the B species 
density. The diffusion rate distributions are rescaled by to the 
same average for better comparison. 

displacement distribution. Figure |2H1 shows the A species 
displacement distributions Ua(x, t) measured at differ- 
ent times from most of the systems examined thus far 
(apart from the case of radioactive B decay), all scaled 
to have approximately the same second moment {x\). 
Aside from the case of one-dimensional B pair annihi- 
lation with its strong anticorrelations, all distributions 
appear to at least approximately collapse to the same 
scaling function. As we have seen previously, compared 
with a Gaussian this common distribution features an 
excess of particles both in its 'tails' and peak. For the 
one-dimensional pair annihilation case these deviations 
are markedly enhanced. We find this apparent univer- 
sality in the A particle displacement distribution quite 
surprising. Yet in each case depicted in Fig. [251 cor- 
relations developed between the B particles, leading to 
significant inhomogeneities in their spatial distribution. 
However, assuming that the A-B coupling is suffi- 



ciently weak that approximate A spatial homogeneity 
is maintained, the distribution of the time-averaged ef- 
fective A diffusivities is constrained by the global value 
D(pB(t)). We compared these quantities for a one- 
dimensional DP system (system size 10000, 10 runs with 
A = 0.01) and indeed saw excellent agreement at long 
times (see the inset of Fig. 1271 , indicating that at each 
time step on average one A particle hops to each B site. 
The discrepancy at small times can be attributed to nu- 
merical errors in computing the integral. Apparently the 
time-averaged diffusion rate distribution over the A par- 
ticles is common to the bulk of the reaction systems stud- 
ied here. To address the apparent universality in systems 
of positive or weak (but existing) B species correlations, 
we should note that localization, by which we mean the 
annihilation of the nearest B particle to the newly lo- 
calized A, is often impermanent when the A resides in 
a B cluster that results from positive correlations. Fig- 
ure [S7| shows the effective A diffusion rate distribution at 
various times for the one-dimensional DP system. The 
distribution at t = 10000 is still fairly peaked, which may 
be accounted for by the initial high density of B parti- 
cles leading to homogenization of the diffusion rates. The 
distributions at later times (scaled to match the average 
value of the t = 10000 distribution) seems approximately 
constant, though perhaps the peak is slowly shifting to- 
wards smaller values. Finally, we note that the measured 
values of D were approximately the same for all systems, 
namely D « 0.5 in d = 1 and D rj 0.25 in d = 2. 
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